Integrated analysis of waveguide dispersed GPR pulses using deterministic and Bayesian inversion methods

Ground-penetrating radar (GPR) data affected by waveguide dispersion are not straightforward to analyse. Therefore, waveguide dispersed common midpoint measurements are typically interpreted using so-called dispersion curves, which describe the phase velocity as a function of frequency. These dispersion curves are typically evaluated with deterministic optimization algorithms that derive the dielectric properties of the subsurface as well as the location and depth of the respective layers. However, these methods do not provide estimates of the uncertainty of the inferred subsurface properties. Here, we applied a formal Bayesian inversion methodology using the recently developed DiffeRential Evolution Adaptive Metropolis DREAM(ZS) algorithm. This Markov Chain Monte Carlo simulation method rapidly estimates the (non-linear) parameter uncertainty and helps treat the measurement error explicitly. We found that the frequency range used in the inversion has an important influence on the posterior parameter estimates, essentially because parameter sensitivity varies with measurement frequency. Moreover, we established that the measurement error associated with the dispersion curve is frequency dependent and that the estimated model parameters become severely biased if this frequency-dependent nature of the measurement error is not properly accounted for. We estimated these frequency-dependent measurement errors together with the model parameters using the DREAM(ZS) algorithm. The posterior distribution of the model parameters derived in this way compared well with inversion results for a reduced frequency bandwidth. This more subjective method is an alternative to reduce the bias introduced by this frequencydependent measurement error. Altogether, our inversion procedure provides an integrated and objective methodology for the analysis of dispersive GPR data and appropriately treats the measurement error and parameter uncertainty. separation is increased while keeping the same midpoint. In such a CMP measurement, reflected GPR waves can be identified by their hyperbolic shape, which can be used to estimate the depth of the reflecting layer and the average propagation velocity of the layer above the reflecting boundary (Greaves et al. 1996; van Overmeeren et al. 1997; Dannowski and Yaramanci 1999; Endres et al. 2000; Bohidar and Hermance 2002; Garambois et al. 2002; Grote et al. 2003; Lunt et al. 2005; Turesson 2006; Gerhards et al. 2008). Another wave that has been used for soil water content determination is the direct subsurface ground wave, which is the direct transmission from the transmitter to the receiver antenna through the top part of the soil. This ground wave can be recognized in a CMP measurement by its traveltime curve that shows a linear increase in arrival time with antenna separation. The ground wave propagation velocity can be determined from the slope of the traveltime curve. The ground wave velocity has been widely used to measure the spatio-temporal development of soil water content variability (e.g., Galagedara et INTRODUCTION Ground-penetrating radar (GPR) is a geophysical technique that uses pulsed electromagnetic waves to explore the subsurface. The transmitted waves will be partly reflected and partly transmitted when contrasts in dielectric permittivity associated with subsurface structures occur. The propagation velocity of the GPR waves depends on the dielectric permittivity, which in turn can be related to soil moisture content and soil porosity amongst other factors (e.g., van Overmeeren et al. 1997; Huisman et al. 2001; Galagedara et al. 2003, Huisman et al. 2003a, Moysey 2004, Bradford 2008, Westermann et al. 2010; Haarder et al. 2011; Rhim 2011; Steelman and Endres 2011). For on-ground GPR, propagation velocity and therewith, the dielectric permittivity can be determined when GPR measurements are made with multiple antenna offsets, for example using a common midpoint (CMP) measurement where the antenna J. Bikowski et al. 2 © 2012 European Association of Geoscientists & Engineers, Near Surface Geophysics, 2012, 10, xxx-xxx error in the model, which is often introduced by simplifying assumptions regarding the modelling of GPR and the representation of the subsurface. Clearly, it is desirable to simultaneously estimate the ‘best’ model parameters and their associated uncertainty. Bayesian inversion algorithms based on Markov Chain Monte Carlo (MCMC) simulation methods are particularly well suited for this task. Such methods are not new in the field of geophysical inversion (Mosegaard and Tarantola 1995; Sambridge and Mosegaard 2002) but the ever increasing computational power and development of advanced MCMC simulation schemes has resulted in their increased use in recent years, especially in the field of hydrogeophysics (e.g., Strobbia and Cassiani 2007; Irving and Singha 2010; Hinnell et al. 2010; Huisman et al. 2010). The aim of this study is to present an integrated analysis of parameter uncertainty associated with the inversion of synthetic and experimental GPR data with waveguide dispersion using MCMC simulation methods. In order to do so, we first describe deterministic and MCMC inversion methods for dispersive GPR data. Then, we use synthetic dispersive GPR data to illustrate the effects of frequency bandwidth and measurement noise on the al. 2003; Huisman et al. 2003b; Weihermüller et al. 2007). Although the ground wave is promising for soil water content measurements, difficulties arise when the subsurface is highly heterogeneous with distinct layers or gradients in soil moisture content that introduce thin horizontal layers with a strong contrast in dielectric permittivity. If the depth of these layers is comparable to or smaller than the wavelength of the GPR signal, they can act as a waveguide in which the electromagnetic waves are trapped. This leads to positive interference related to total reflection of the trapped wave at the boundaries of the layer. Field situations where such waveguides have been reported include a thin ice sheet floating on water (Arcone 1984; van der Kruk et al. 2007), an organic-rich sandy silt layer overlying a gravel unit (Arcone et al. 2003; van der Kruk et al. 2006), a mountain slope with a soil cover (Strobbia and Cassiani 2007) and thawing of a frozen soil layer (van der Kruk et al. 2009; Steelman et al. 2010). In the presence of such waveguides, CMP measurements are difficult to interpret because the arrival time and the first cycle amplitude of the ground wave cannot be identified due to interfering waves and dispersion. To enable interpretation of GPR data with waveguide dispersion, van der Kruk et al. (2006) presented a deterministic inversion algorithm to simultaneously estimate the thickness and permittivity of the dispersive waveguide and the permittivity of the material below the waveguide. This work was inspired by inversion algorithms for Rayleigh and Love waves commonly observed in multi-offset seismic data and were extended to higher order modes by van der Kruk (2006, 2007). More recent extensions include the inversion for multiple layers acting as waveguides (van der Kruk et al. 2010). Inversion of GPR data affected by waveguide dispersion requires a forward model that accurately describes the dispersive characteristics of GPR data for a given subsurface structure described by a set of model parameters. Optimization methods are then used to seek a set of model parameters that minimizes the discrepancy between simulated and measured GPR data. In general, numerical modelling and inversion methods for GPR data have greatly improved in the last decade, which obviously enhances the inversion quality, as well as the range of GPR data that can be inverted. Yet, traditional GPR inversion algorithms are deterministic and estimate only a single ‘best’ set of model parameters without consideration of parameter uncertainty (e.g., Pettinelli et al. 2007; Steelman and Endres 2010; Wollschläger et al. 2010). Therefore, it is not yet well established how errors in GPR measurements and models propagate through the processing and inversion of dispersive GPR data. For non-dispersive data, reported confidence intervals of wave velocity and hence, implicitly of dielectric permittivity, vary widely depending on the field settings and methods used (e.g., Jacob and Hermance 2004). Typical errors of GPR measurements introduced during data acquisition include inaccuracies in offset, timing, antenna orientation and other antenna effects (Slob 2010). An additional error source that is more complicated to address is the structural FIGURE 1 Schematic outline of the analysis approach. The measured CMP data are transformed into a dispersion curve that is affected by the measurement error. The forward model provides the modelled dispersion curve given the model parameters describing the waveguide. Either a deterministic or Bayesian inversion method is applied. Only the Bayesian inversion approach can be used to simultaneously estimate model parameters and the measurement error associated with the dispersion curve. The result of the deterministic method is a single value for each model parameter, whereas the Bayesian inversion with MCMC simulation (DREAM(ZS)) returns the posterior probability distribution of the model parameters, which provides the most likely model parameters and their uncertainty. Integrated analysis of waveguide dispersed GPR pulses 3 © 2012 European Association of Geoscientists & Engineers, Near Surface Geophysics, 2012, 10, xxx-xxx velocity in the measured dispersion curve. This also suggests that an adequate representation of low frequencies is required to obtain accurate estimates of ε2. For high frequencies, the dispersion curve asymptotically approaches , which suggests that an adequate representation of high frequencies is required for accurate estimates of ε1. However, the direct estimation of ε1 and ε2 using the simple procedure outlined above can be difficult because of the low signal-to-noise ratio in the dispersion curve for low and high frequencies. DETERMINISTIC INVERSION Two inversion algorithms are used in this study: the deterministic inversion algorithm of van der Kruk et al. (2006) and the DiffeRential Evolution Adaptive Metropolis, DREAM(zs) algorithm of Vrugt et al. (2009). The deterministic inversion approach of van der Kruk (2006) aims to find the model parameters (m = ε1, ε2, h) that minimize the mean of the absolute difference between modelled and measured dispersion curves. In contrast to van der Kruk et al. (2006), we describe the misfit, C, between the measured dispersion curve, vmeas(f) and modelled dispersion curve, vmod(f,m), using a mean squared difference,